Creep modeling of composite materials based on improved gene expression programming

In this article, a new method for creep modeling and performance prediction of composite materials is presented. Since Findley power-law model is usually suitable for studying one-dimensional time-dependent creep of materials under low stress, an intelligent computing method is utilized to derive three temperature-related sub-functions, the creep model as a function of time and temperature is established. In order to accelerate convergence rate and improve solution accuracy, an improved gene expression programming (IGEP) algorithm is proposed by adopting the probability-based population initialization and semi-elite roulette selection strategy. Based on short-term creep data at seven temperatures, a bivariate creep model with certain physical significance is developed. At fixed temperature, the univariate creep model is acquired. R2, RMSE, MAE, RRSE statistical metrics are used to verify the validity of the developed model by comparison with viscoelastic models. Shift factor is solved by Arrhenius equation. The creep master curve is derived from time–temperature superposition model, and evaluated by Burgers, Findley and HKK models. R-square of IGEP model is above 0.98 that is better than classical models. Moreover, the model is utilized to predict creep values at t = 1000 h. Compared with experimental values, the relative errors are within 5.2%. The results show that the improved algorithm can establish effective models that accurately predict the long-term creep performance of composites.

www.nature.com/scientificreports/ Findley model are derived from improved gene expression programming algorithm to establish bivariate creep model. Compared with classical viscoelastic models, the validity of univariate model is verified by four statistical metrics at fixed temperature. Lastly, creep master curve is drawn from time-temperature superposition model based on shift factor. The developed model is applied to predict long-term creep performance of composites so that high prediction accuracy of the model is validated.

The proposed methodology
In this section, the flow chart of overall research methodology is given as shown in Fig. 1. The methodology is divided into two stages: creep modeling and performance prediction. Further detail steps of the research are discussed in the subsequent subsections.
Preparation of specimen. The matrix material of this experiment is m-benzene type unsaturated polyester resin FC518, which was supplied by Shanghai Fuchen Chemical Co., Ltd. The reinforcement materials are made up of alkali-free glass fibers, with the specifications of winding yarn 2400 Tex and chopped strand mat 450 g/m2, which were provided by Hebei Zhongyi Composite Materials Co., Ltd. The experimental specimens are: resin (R), fiber chopped strand mat (CSM) and fiber circumferential winding (FWC). According to the standard GB/T 1449-2005, INSTRON5828 is used to test the initial flexural strength (σ) of specimens. The resin mass content (W) of each specimen is tested based on the standard GB/T 2577-2005, and the results are given in Table 1. The size of each specimen is determined by the above-mentioned standard, the thickness h = 5 mm, the width b = (2.5 ± 0.5) h, and the length L = (18 ± 2) h. The constant load applied by INSTRON5848 universal testing machine is 20% of the initial flexural strength, and these testing data are automatically read by computer with a time interval of 0.1 s.   www.nature.com/scientificreports/ Overview of gene expression programming. Gene Expression Programming (GEP) invented by Ferreira is derived and improved from genetic algorithm and genetic programming, it is an efficient tool for developing models and consists of chromosome with fixed length. Each gene in the chromosome contains a head h and a tail t , there exists the following relationship:t = h(n − 1) + 1 , n is the total number of arguments within a function(maximum arity). The head of each gene contains both function symbols and terminal symbols (e.g.{+ , −, *, /,√,cos, tan, log, 6, x, a, b}). While the tail only contains terminal symbols that are made up of constants and variables (e.g.{8, y, c, d}). The chromosomes can be viewed as genomes that are modified through selection, crossover, mutation, transposition and recombination operations. GEP is developed based on two essential elements: chromosome and expression tree (ET). The genotype of GEP is chromosome, and the phenotype is ET that is composed of nonlinear entities with different sizes and shapes. For example, the chromosome consists of one gene, and the genotype of individual is: * − sinQ + cab/bababbaaba, the part in bold is the tail. The gene has a head length of 9 and a tail length of 10, so the total length of gene is 19. The genome and expression tree can be converted into each other in a certain way, as shown in Fig. 2.
The mathematical equation corresponding to genotype can be expressed as: √ a− b + b a * (sinc) . Simultaneously, the fitness value Fitness(i) of an individual i is calculated, as given in the Eq. (1).
where M is the selected range,C(i, j) is the value returned by an individual i for fitness case j(out of n fitness cases),T(j) is the target value for fitness case j . If C(i, j) = T(j) , there is Fitness(i) = n · M , the system can find the optimal model for itself by this way. Therefore. GEP greatly surpasses existing adaptive techniques 37 .

Proposition of improved GEP.
The individuals of GEP have linear genotype and non-linear phenotype. Simultaneously, GEP not only contains all the genetic operators of traditional evolutionary algorithm but also introduces some new operators, which brings some challenges to the convergence rate and solution accuracy. Although GEP algorithm has flexible encoding/decoding methods and evolutionary operations, when there are many terminal symbols in the head of gene, it is easy to generate invalid individuals; when the fitness function is selected, the lack of population diversity results in slow convergence, and it is easy to fall into local optimum. Therefore, this article proposes an improved GEP (IGEP) algorithm. The individuals are initialized by probability to accelerate convergence rate; the semi-elite roulette selection is performed to improve solution accuracy. Its flow chart is shown in Fig. 3.
The detailed steps of the algorithm are given in Table 2.

Algorithm complexity analysis.
The low time cost of IGEP algorithm is very important to be used to build the model. Given the maximum number of iterations is MAXGEN , the size of population is N , the size of elite population is M , the length of gene is len , and the size of sample data is S . As can be known from the algorithm, in Step 1, individuals with a length of len are traversed and gene encoding is performed. Therefore, the time complexity of population initialization process is O(N · len) . In Step 2, the fitness value of each individual is evaluated, so the time complexity is O(N · S) . In Step3, firstly, the ratio of individual fitness to overall fitness is calculated, and its time complexity is O(N) ; secondly, the semi-elite roulette strategy is utilized to select individuals, the time complexity is O N 2 ; thirdly, the sorting algorithm is employed to select elite population, and its time complexity  www.nature.com/scientificreports/ genes are exchanged, its time complexity is O(N · len) . In summary, the time complexity required for one iteration is O(3N · len + N · S + N 2 + N + N · log(N)) . After removing constant term and simplifying the formula, the total time complexity of all iterations is O((len + S + N + log(N)) · N · MAXGEN) 38 .
The physical model. Burgers model. Burgers model is a combination of Maxwell and Kelvin-Voigt elements, it is one of the most widely used models to give the relationship between morphology of composites and their creep behavior 39 , which is a four-element model, as shown in Fig. 4.
For the most general case of linear viscoelastic materials, the total creep strain is essentially the sum of three separate parts: ε 1 is the instantaneous elastic deformation; ε 2 is the delayed elastic deformation; ε 3 is the Newtonian flow, it is the same as the deformation of a viscous liquid that obeys Newton's law of viscosity. The total strain ε B (t) as a function of time corresponds to the following Eq.  where t is the time,σ 0 is the applied stress, ε H (t) is the total strain, C H (t) is the creep compliance, E i and η j are the model parameters, i=0, 1, 2 , j=1, 2.E 0 is the initial elastic modulus; E 1 and E 2 are the elastic moduli of Kelvin springs, respectively; η 1 and η 2 are the viscosities of Kelvin dashpots, respectively.

The phenomenological model. Findley power-law model. The phenomenological model developed by
Findley introduces a mathematical expression to describe creep behavior of composite materials that is more suitable for the prediction of creep deformation, it can effectively predict mechanical performance of composites. In this model, the creep response can be divided into time-independent and time-dependent strains, creep strain can be expressed as follows: where ε 0 is the initial stress-dependent and time-independent elastic strain,ε c is a coefficient related to stress and temperature, t is the time, n is a stress-independent and temperature-dependent material constant 40 . Under constant stress, the subsequent form (7) could be derived, where C 0 is the initial temperature-dependent creep, m is a temperature-related coefficient, and n is a dimensionless material parameter that is dependent of temperature.
Since the specific mathematical form of Findley model with time and temperature has not been deduced in the theoretical analysis, at different temperatures, C can be determined as a bivariate function of both time and temperature. Therefore, Findley model is considered as a modeling framework, the modified model is represented as: Time-temperature superposition model. Assuming that creep compliance is a function related to time and temperature, the creep behavior of composites at low temperature for a long time can be predicted by using shortterm creep data at high temperatures. The creep compliance curve C T ref , t/φ T at reference temperature T ref can be constructed by shifting short-term compliance curve C(T i , t) at different temperatures along the logarith- Figure 4. Schematic diagram of Burgers model. www.nature.com/scientificreports/ mic time axis by shift factor φ T , and so the smooth creep master curve is derived, which is time-temperature superposition (TTSP), the calculation equation is as follows: Assuming the activation energy is constant, time-temperature shift factor φ T is obtained to construct creep master curve, it is in good quantitative agreement with the Arrhenius equation, the formula is given in (9), which provides a reliable method for predicting long-term creep performance of composite materials.
where E a is activation energy [ kJmol −1 ], R is the universal gas constant with a value of 8.314 × 10 −3 kJK −1 mol −1 ,T is the testing temperature [K]. Equation (9) is applicable for temperature below glass transition temperature.

Creep data and experimental settings
Data description. Three-point bending tests are carried out under constant load. The temperatures of R, CSM and FWC specimens are set to 20 °C, 25 °C, 30 °C, 35 °C, 40 °C, 45 °C, and 50 °C, respectively. According to the standards, these specimens need to be maintained within a constant temperature chamber for 20 min before testing to ensure that the experimental temperature is reached. Short-term (1 h) flexural creep performances of three specimens are tested at seven temperatures, and creep data ranging from 0 to 3600 s are obtained. The resin content of R specimen is 100%, without any constraint of reinforcement materials, so the creep compliance and creep growth rate are the largest, and its creep resistance is the weakest; the resin content of FWC specimen is the lowest, and its continuous fibers have the strongest constraint effect on resin deformation, so its creep compliance is the smallest and creep resistance is the strongest; the resin content of CSM specimen is relatively high, many interfaces lead to stress concentration, and the constraint effect of chopped fibers on resin is not as strong as that of continuous fibers, so its creep compliance and creep resistance are between the two. Therefore, the creep compliance C-time t curves of R, CSM and FWC specimens could be drawn, as shown in Fig. 6.
Experimental settings. Various parameters are involved in the establishment of IGEP model, and affect generalization capability of the model. In order to get a more accurate IGEP model and reduce the time complexity, appropriate parameters need to be set for problem solving, including fitness function, the number of iterations, population size, the number of genes, linking function and probabilities of genetic operators. Based on multiple trials, the final parameters selected for IGEP algorithm are given in Table 3.
Evaluation metric. Four evaluation metrics, namely, coefficient of determination R-squared R 2 , root mean square error RMSE, mean absolute error MAE and relative square root error RRSE, are used to evaluate the performance and compare prediction accuracy of models. These criteria are calculated as follows: Figure 6. Creep compliance C-time t curves of three specimens at seven temperatures. www.nature.com/scientificreports/ where n is the number of data points, y i is the measured value, y is the average value, and ŷ i is the predicted value. R 2 measures the degree of correlation, the larger the value of R 2 , the better the performance of model; RMSE is a measure of the residual variance, lower RMSE represents more accurate estimation; the smaller the values of MAE and RRSE, the better the performance of model.

Creep modeling results and validation of the model
The experiments are implemented on a PC with Intel Core i5-4460 3.20 GHz CPU, 8 GB memory, Win7 64-bit operating system, and the software environment is MATLAB R2016a.
Time-temperature bivariate creep modeling results. Based on short-term creep data from 0 to 3600 s, Findley model is modified to be expressed as a function of time and temperature by IGEP algorithm. Therefore, the time-temperature bivariate IGEP models for three specimens are established, the modeling results are given in Table 4, where a i (i = 1, 2, 3,…) is the model parameter, C 0 , m and n are three sub-functions related to temperature T respectively, and R 2 of three models are above 0.  Table 4. IGEP creep models for three specimens.      Table 8. Similarly, IGEP models for CSM and FWC specimens are analyzed by dimension reduction, and the creep curves at 50 °C are obtained. Compared with viscoelastic models, the curve fitting results are plotted in Fig. 9b,c. Simultaneously, the overall performance of IGEP model can be validated by four metrics R 2 , RMSE, MAE and RRSE, the values are provided in Table 8.
Since most of creep models are time-related univariate models, and there are few models with multiple variables, a new bivariate modeling program is developed by IGEP in this work, the effect of temperature is introduced into the traditional Findley power-law creep equation. It can be clearly seen from the table that R 2 values of univariate IGEP model for three specimens are above 0.92 by dimension reduction analysis, the coefficient of determination of four models are relatively high and close to each other. The results show that the fitting curve of IGEP model is almost in good agreement with experimental data.

Time-temperature superposition creep modeling results. Calculating activation energy is a very
useful technique to estimate shift factor for time-temperature superposition without constructing complete master curve. The activation energies E a of R, CSM and FWC specimens are obtained by dynamic mechanical thermal analysis to be 365.50 kJ/mol, 337.07 kJ/mol and 319.66 kJ/mol, respectively. Assuming E a is valid only below material's glass transition temperature. In this article, 23 °C is selected as reference temperature T ref . Since some experimental temperatures are higher than reference temperature 23 °C, others are lower than 23 °C. For T > T ref , the logarithm of shift factor lg φ T is negative resulting in right-shifted creep compliance www.nature.com/scientificreports/ curve. On the contrary, for T < T ref , the logarithm of shift factor lg φ T is positive resulting in left-shifted creep curve. According to Arrhenius equation, the logarithm of shift factor for three specimens are calculated as given in Table 9. It is clearly seen that the order of lg φ T for three specimens at the same temperature is as follows: lg(R) > lg(CSM) > lg(FWC) . The larger the logarithm of shift factor, the greater the effect of temperatures on creep performance of composites. Therefore, the sensitivity of creep to temperatures for three specimens is: R > CSM > FWC. When short-term experimental data of creep compliance C-time t at seven temperatures are used, creep master curve of R, CSM and FWC specimens can be derived from TTSP, as shown in Fig. 10. Findley model is a parametric phenomenological model suitable for creep behaviour under low stress conditions. Burgers model and HKK model are classical physical models. At present, there are different methods for master curve fitting. The abscissa axis in Fig. 10 represents the logarithm of time lg t . In order to facilitate the observation, the abscissa axis is converted into time t, and is plotted to provide reference for engineering structural design. IGEP model and viscoelastic models are established by fitting the data on creep master curve for three specimens. The results are shown in Fig. 11. Moreover, the metric values of four models are calculated, as given in Table 10.
Apparently, it can be revealed that R 2 of IGEP model is above 0.98, and the curve fits well with experimental data. R 2 , RMSE, MAE and RRSE are used as evaluation metrics, the fitting effect of IGEP model is better than that of Findley model, and far better than that of Burgers model and HKK model, indicating that IGEP model can well describe long-term creep performance of composite materials. When the creep compliance values of R and CSM specimens are enlarged by 10 10 times, and creep compliance values of FWC specimen are enlarged by 10 11 times, the model parameters are figured out with the help of computational software Origin 2018, the results are provided in Tables 11, 12 and 13.  Figure 10. Creep master curves for three specimens. www.nature.com/scientificreports/ Figure 11. Creep master curves and fitting curves for three specimens.  When the creep compliance at t = 1000 h is selected for analysis, four models are established by fitting master curve to predict the values at t = 1000 h. The predicted value of TTSP model is compared with experimental value at 23 °C, and then the relative error δ TTSP is calculated, the results are provided in Tables 14 and 15. It can be seen that relative error δ TTSP predicted by TTSP model for R specimen is 5.18%; relative error δ TTSP for CSM specimen is 2.22%; and relative error δ TTSP for FWC specimen is 1.15%, all are within 6%. It is well proven that long-term flexural creep life of composites can be accurately predicted through an accelerated testing method at high temperatures.
However, it is clearly seen from Table 14 that the prediction effect of IGEP model for R specimen is better than that of TTSP model and Findley model, far better than that of Burgers model and HKK model; the prediction effect of IGEP model for CSM specimen is comparable to that of Findley model, better than that of TTSP model, and far better than that of Burgers model and HKK model; the prediction effect of IGEP model for FWC specimen is better than that of TTSP model, and comparable to that of other creep models. It is concluded that IGEP model is a better way to simulate creep master curve.
At the same time, taking relative error δ between creep value predicted by each model and experimental value at t = 1000 h as a statistical metric, it can be seen from Table 15 that the relative error δ IGEP of IGEP model for R specimen is the smallest, it is 5.11%; the relative error δ IGEP of IGEP model for CSM specimen is almost the same as the error δ Findley of Findley model, it is 0.61%; at t = 1000 h, the prediction effect of each model for FWC specimen is better than that of TTSP model, and the relative error δ is very small, all are below 0.6%. The predicted values are extremely consistent with experimental values. The experiments and theory are integrated to verify the validity of accelerated characterization method. The comparison of developed models and accelerated testing results indicates that IGEP model has better prediction accuracy than Burgers, Findley and HKK models in describing long-term creep performance of composite materials.
Discussion. Creep modeling of composite materials is a subject widely studied in the field of material science and engineering. This article only investigates the effect of time and temperature on flexural creep behavior of composites that is very important to the service life. However, under the complex conditions, there are many factors involved in creep failure of materials, such as humidity, atomic migration and diffusion, crack initiation and propagation, fiber morphology and orientation, there exists uncertainty in creep properties of composites, so that the empirical prediction model is not accurate enough. In addition, the joint effect of various factors makes it difficult to simulate the evolution process of creep from a microscopic perspective. Therefore, a swarm intelligent algorithm can be utilized to establish mathematical relationship model between multiple factors and output from a macroscopic perspective. The randomness and fuzziness of creep are not considered that results in the failure of classical models. So the fuzzy random method is used to improve the traditional particle swarm algorithm to obtain an efficient model to describe creep performance 41 . The operation of instrument and the test of specimen result in certain errors in the data obtained. The creep test under constant load is performed, but the load is variable in practical application. An effective model for describing creep properties of composites under step loading and unloading conditions is established to provide theoretical support for deformation analysis and long-term stability 42 .
The intelligent evolutionary algorithm is easy to implement and has strong scalability by selecting different basis functions such as exponential function and power function. When there are few experimental samples, the useful information can still be analyzed and extracted from the data, so that the testing workload in the process www.nature.com/scientificreports/ of creep modeling is reduced. The modeling of alternative methods proves that the prediction of machine learning algorithm is superior to other methods in the literature 43 , it has wider engineering applicability and higher prediction accuracy in describing long-term creep performance of composites. GEP is an efficient evolutionary algorithm, it can be regarded as a promising approach to devise empirical models based on experimental phenomena and variation laws, Since the creep experiments at room temperature need to take a long time, the application of accelerated characterization method can reduce its time cost by short-term creep data. Although mechanical testing is one of the most direct ways to study mechanical properties of materials, the time-consuming and sophisticated creep tests could be avoided through computer simulation using GEP.

Conclusions
To summarize this article, an intelligent computing method is proposed for creep modeling of composite materials. In order to derive three temperature-related sub-functions of Findley model, an improved GEP algorithm is developed to establish bivariate model. The probability-based population initialization and semi-elite roulette selection are adopted to accelerate convergence rate and improve solution accuracy. Moreover, compared with Burgers, Findley and HKK models, the validity of univariate model at fixed temperature is verified by R 2 , RMSE, MAE and RRSE metrics. Lastly, the short-term creep curves are plotted as creep master curve based on shift factor, the relative error at t = 1000 h is used as a statistical metric. IGEP model established by fitting master curve has lower prediction errors for three specimens, all are within 6%. The experimental results indicate that IGEP model can accurately predict long-term creep performance of composite materials. This work not only expands the application field of GEP algorithm, but also provides a new method for creep modeling.
In future work, except for TTSP, other superposition models could be extended and are reasonably studied to accelerate the characterization of long-term performance. When the effect of fiber content and surface treatment on creep properties of composites is studied further, GEP algorithm would be efficiently utilized to develop multivariable creep model as a function of temperature, stress and fiber, which is of great significance to investigate creep behavior for the design and life prediction.

Data availability
The data used to support the findings of this study are available from the corresponding author upon reasonable request.